*------------------------------------------------------------------------------*
			*** FIGURE S11a. POLLUTION IN MUNICH 2016 SUNDAYS ***
*------------------------------------------------------------------------------*

use "${rawdata}pollution/pollution_UBA_d_19902019.dta", clear

*------------------------------------------------------------------------------*

*Keep year 2016
keep if year == 2016

*Keep relevant variables
keep station year month day PM10_ug

*Gen date and day of the week
gen date = mdy(month, day, year)
format date %td
gen dow = dow(mdy( month, day, year))

*Keep only Sundays
keep if dow == 0

*Keep Bavarian stations
keep if substr(station, 1, 4) == "DEBY"

*Add station coordinated
merge m:1 station using "${rawdata}pollution/meta_station_UBA.dta", keepusing(latitude longitude) keep(master match) nogen

*Add Munich coordinates
gen latitude_kreis = 48.13642
gen longitude_kreis = 11.57755

*------------------------------------------------------------------------------*

**Distance
distance latitude_kreis longitude_kreis latitude_d longitude_d , a(kreis) b(station)

*Keep stations within 30km
keep if dist_KRST <= 30

*Keep stations with valid PM10
keep if PM10_ug != .

**Inverse-distance weighted PM10
*Create weights
gen double weight = 1/dist_KRST
bys date: egen double weight_tot = total(weight)
replace weight = weight/weight_tot
drop weight_tot

*Weighted mean
bys date: egen double PM10_weight = total(PM10_ug*weight)
bys date: keep if _n==1

rename PM10_weight PM10_weight_30km

*Gen PM10 in tens of micrograms
gen PM10_weight_10ug = PM10_weight_30km / 10

*------------------------------------------------------------------------------*

**Graph pollution 2017 sundays
drop if date == td(1jan2017)  // outlier
graph twoway (line PM10_weight_10ug date), ytitle(PM10 (10 {&mu}g/m{sup:3})) xtitle("") ///
		tlabel(15jan2016 "Jan" 15feb2016 "Feb" 15mar2016 "Mar" 15apr2016 "Apr" 15may2016 "May" ///
		15jun2016 "Jun" 15jul2016 "Jul" 15aug2016 "Aug" 15sep2016 "Sep" 15oct2016 "Oct" 15nov2016 "Nov" 15dec2016 "Dec") ///
		graphregion(fcolor(white)) name(sundays, replace)

graph export "${outputs}figure_S11a.pdf", as(pdf) replace

*------------------------------------------------------------------------------*

clear

exit
